%% The upconingbox
% - wells sheet


clear variables; close all; clc;

basename    = 'Upconing';
save basename basename;

%% Setup the grid
ML  = 2250;	% [m] Length of model;    x-ax
MW  = 40;   % [m] Width of the model; y-ax
MDtop = -25;% [m] Depth of the model; z-ax
MDbot = -150;% [m] Depth of the model; z-ax
dx  = 50;	% [m] grid cell width
dy  = 20;	% [m] grid cell length = thickness of model; The model has become a 2D model, Hele Shaw
dz  = 10;	% [m] grid cell size vertical
depthSaltFront = -100; % [m] The depth of the salt front



xGr = [0:dx:ML ML];
yGr = [0:dy:MW MW];
zGr = [MDtop:-dz:MDbot MDbot];

gr  = gridObj(xGr, yGr, zGr);

%% Get period data
[perdata.pernams,perdata.pervals]   = getPeriods(basename);
[perdata.perlen]            = perdata.pervals(:,strmatchi('PERLEN', perdata.pernams));
[perdata.nper]              = length(perdata.perlen);
[perdata.time]              = cumsum(perdata.perlen);

%% Set initial variables
freshwater  = 0; % relative density compared to saltwater defined in SEAWAT-sheet
saltwater   = 1; % relative density compared to saltwater defined in SEAWAT-sheet

IBOUND  = gr.const(2);  % Define active cells

K       = 25;   % [m/day] conductivity
HK      = gr.const(K);  % [m/day] horizontal conductivity
VK      = HK;           % [m/day] vertical conductivity

p       = 0.38;         % [] effective porosity
PEFF    = gr.const(p);

STRTHD  = gr.const(-25);  % [m] is top of the model

ICBUND  = gr.const(saltwater);
STCONC  = gr.const(saltwater);

% add freswater start concentration
iDepthSaltFront = round(interp1(gr.zm,1:length(gr.zm),depthSaltFront));
STCONC(:,:,1:iDepthSaltFront) = freshwater;

%% Define cellproperties/apply packages
CHDzone = 3;                % needed for gr.bcnZone
IBOUND(:,end,:) = CHDzone;  % apply chd package

% done trough the per-sheet
RECHzone= 6;
IBOUND(:,:,1) = RECHzone; % RECH package

%% Constand head package
CHDDENSOPT      = 2;        % env. hd at ocean boundary, Langevin et al 2008, p22
[CHD, PNTSRC]   = gr.bcnZone(basename,'CHD',IBOUND,{CHDzone 'CHD' 'CHD' CHDDENSOPT},saltwater);

      
%% Well package
[well,WEL,PNTSRC(end+1:end+perdata.nper),NPER]=mf_setwells(basename,gr,HK,'wells');

%% Rech package
% done trough the per-sheet
% simple since rech is constant over time
% RECHintensity = 0.05;    % [m/day]
% % RECH    = gr.const(0);
% % RECH(IBOUNDtop==6) = RECHintensity;
% RECH_ = (squeeze(IBOUND(:,:,1))==6).*RECHintensity;  % Only upper layer in 2d for rech
% RECH=repmat(RECH_,[1, 1, perdata.nper]);

save underneath well gr perdata